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ABSTRACT  OF  THESIS 


THE  INVERTIBILITY  PRINCIPLE  FOR  A  SIMPLE  HURRICANE  MODEL 

\ 

^  i  f  *<''■£- 

We  present  a  reinterpretation  of  Ooyama’s  classic  tropical  cyclone  model  ^in  terms 

of  the  more  recent  theoretical  notions  of  the  potential  thickness  equation,  the  potential 
radius  coordinate  and  the  invertibility  principle.  This  helps  place  tropical  cyclone  theory 
in  a  theoretical  framework  closer  to  that  of  midlatitude  theory.  We  first  present  a  shallow 
water  model  of  axisymmetric,  frictionless  flow  of  homogeneous  incompressible  fluid  on  an 
/-  plane.  Then-we4ntroduc^  the  potential  thickness,  the  inverse  of  potential  vorticity, 
c(  and  the  equation  for  its  evolution.  We  transform  the  system  from  physical  space 
to  absolute  angular  momentum  or  potential  radius  space.  This  eliminates  the  radial 
component  of  the  wind  from  the  problem  and  provides  better  resolution  in  areas  of  large 
vorticity.  We  derive  and  solve  the  invertibility  principle  using  five  different  methods- 
solving  for  a  potential  function  using  the  shooting  method,  the  fluid  depth  using  the  same 
method,  the  fluid  depth  using  a  nonlinear  equation  solving  “black  box”,  a  transformed 
velocity  using  a  tridiagonal  matrix  equation  solver,  and  the  transformed  velocity  with  the 
nonlinear  equation  solver.  The  first  four  methods  can  not  be  generalized  to  two  layers.  A 
short  review  of  Ooyama’s  model  shows  how  the  wind  field  is  determined  from  the  radial 
mass  flux.  Then  we  generalize  the  concepts  of  the  shallow  water  case  to  a  two  layer  model 
using  the  nonlinear  equation  solver  to  solve  the  invertibility  for  the  transformed  velocity. 

Robert  M.  Fogarty 

Department  of  Atmospheric  Science 

Colorado  State  University 

Fort  Collins,  Colorado  80523 

Summer  1989 
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Chapter  1 


INTRODUCTION 

When  studying  atmospheric  dynamics  it  is  often  useful  to  begin  with  the  simplest 
form.  By  reducing  a  problem  to  its  least  complicated  terms  we  can  learn  about  its  most 
important  features.  In  1969  Ooyama  presented  what  is  apparently  the  simplest  closed  set 
of  equations  that  describe  the  life  cycle  of  tropical  cyclones.  His  model  has  the  simplest 
geometry,  moist  physics,  and  the  minimum  vertical  resolution.  Ooyama’s  model  can  even 
be  thought  of  as  the  tropical  meteorologist’s  analogue  of  the  mid-latitude  meteorologist’s 
Eady  model  which  is  the  simplest  model  of  baroclinic  instability  (Haltiner  and  Williams, 
1980). 

The  vortex  is  axisymmetric  and  described  on  an  /-plane  (/,  the  Coriolis  parameter,  is 
constant).  The  vertical  structure  is  a  shallow  diagnostic  boundary  layer  and  two  prognostic 
main  layers  above.  Moisture  is  carried  only  in  the  boundary  layer  and  only  deep  convective 
clouds  are  parameterized.  The  simplicity  of  the  model  makes  it  attractive  for  use  in 
trying  different  formulations  of  the  physical  processes  involved  in  the  formation  of  tropical 
cyclones. 

Some  recent  advances  in  atmospheric  dynamics  can  be  very  useful  in  the  study  of 
tropical  cyclones.  The  use  of  potential  vorticity  reasoning  (Hoskins  et  al.,  1985;  Schubert 
and  Alworth,  1987)  in  combination  with  absolute  angular  momentum,  or  potential  radius, 
coordinates  (Shutts  and  Thorpe,  1978;  Schubert  and  Hack,  1983)  has  proven  advantageous 
in  this  endeavor.  The  goal  here  is  to  introduce  these  ideas  in  the  framework  of  Ooyama’s 
model.  From  the  theoretical  standpoint  this  work  is  concerned  with  the  solution  of  the 
invertibility  principle. 

A  historical  review  of  the  study  of  potential  vorticity  is  given  by  Hoskins  et  al.  (1985). 
Some  of  the  important  early  work  was  done  by  Rossby.  He  found  (Rossby,  1939)  that  only 
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the  vertical  component  of  vorticity  is  important  for  describing  the  large  scale  atmospheric 
flow.  In  addition,  he  was  able  to  produce  the  barotropic  model  by  assuming  the  absolute 
vorticity  is  conserved  following  the  motion  of  an  air  parcel.  It  is  important  to  note  that 
this  model  describes  two-dimensional  horizontal  motion.  Rossby’s  later  work  (1940)  gave 
us  perhaps  the  simplest  statement  of  potential  vorticity.  He  showed  that  if  h  is  the  fluid 
depth  then  £ /h  —  constant  following  the  motion,  where  £  is  the  absolute  vorticity.  This 
equation  relates  the  creation  of  vorticity  by  stretching/shrinking  of  vortex  tubes  and  the 
horizontal  advection  of  absolute  vorticity  which  are  the  two  processes  that  are  generally 
most  important  in  the  vorticity  budget.  Ertel  (1942),  working  independently,  generalized 
Rossby’s  result  to  three-dimensional,  nonhydrostatic  flow.  He  showed  that  following  the 
motion  of  an  air  parcel,  the  potential  vorticity,  P  =  p-1£  •  V0  =  constant,  where  p  is 
the  density  of  the  fluid  and  6  is  the  potential  temperature.  Thus,  we  have  the  concept  of 
potential  vorticity  which  Hoskins  et  al.  describe  as  “a  potential  for  creating  vorticity  bv 
changing  latitude  and  by  adiabatically  changing  the  separation  of  isentropic  layers.” 

Because  it  is  conserved,  the  potential  vorticity  was  quickly  recognized  as  a  Lagrangian 
tracer.  An  air  parcel  can  be  identified  once  and  for  all  by  its  potential  vorticity.  However, 
this  quantity  is  much  more  useful  than  just  a  tracer.  Hoskins  et  al.  credit  Kleinschmidt 
(1950a, b;  1951,  1955,  1957)  with  the  earliest  discovery  of  this.  He  recognized  that  the 
complete  flow  structure  can  be  acsertained  from  the  potential  vorticity  by  the  inversion 
of  a  Laplacian  type  operator.  This  idea  of  recovering  the  wind  field  from  the  potential 
vorticity  field  is  known  as  the  invertibility  principle.  For  example,  in  the  barotropic 
case  once  the  vorticity  is  known,  the  streamfunction  can  be  calculated,  and  the  winds 
determined.  However,  it  is  not  as  simple  as  it  first  appears.  There  are  three  necessary 
conditions  to  make  invertibility  possible.  The  problem  must  be  solved  globally,  in  reference 
to  a  basic  state,  and  under  some  balance  condition  (geostrophic,  gradient,  etc.).  The 
necessity  of  a  balance  condition  is  probably  the  reason  potential  vorticity  has  not  been 
used  extensively  until  recently  (Schubert  et  al.,  1986). 

Before  looking  at  some  recent  work  using  potential  vorticity  reasoning  to  describe  the 
structure  of  tropical  cyclones  it  is  useful  to  examine  the  potential  radius  coordinate.  In 


1978  Shutts  and  Thorpe  used  the  absolute  angular  momentum  as  the  radial  coordinate  in 
a  gradient  balance  vortex  model.  They  described  it  as  analogous  to  the  use  of  geostrophic 
momentum  coordinates  in  the  semi-geostrophic  equations.  This  formulation  allowed  dy¬ 
namical  effects  associated  with  departure  from  gradient  balance  to  be  included  implicitly 
in  the  same  way  Hoskins  and  Bretherton  (1972)  showed  ageostrophic  motion  is  included 
implicitly  in  the  semigeostrophic  theory  of  frontogenesis.  Schubert  and  Hack  (1983)  trans¬ 
formed  the  Eliassen  (1952)  balanced  vortex  model  to  potential  radius  coordinates  and  used 
it  to  study  tropical  cyclones.  This  model  describes  a  slow,  quasi-balanced,  thermally  or 
frictionally  driven  axisymmetric  vortex.  If  friction  is  assumed  to  be  confined  to  a  shallow 
boundary  layer  then  absolute  angular  momentum  is  conserved  above  it.  This  allows  the 
use  of  potential  radius  above  the  boundary  layer.  They  introduced  this  coordinate  as  a 
dependent  variable  in  Eliassen ’s  equations  then  switched  the  roles  of  actual  radius  and 
potential  radius  making  it  the  independent  variable.  The  result  of  this  transformation 
is  the  stretching  of  regions  of  positive  relative  vorticity  and  the  shrinking  of  regions  of 
negative  relative  vorticity.  This  again  is  compared  to  Hoskins  and  Bretherton ’s  results  in 
their  frontogenesis  study.  These  studies  show  that  the  potential  radius  is  indeed  a  useful 
tool. 

The  question  arises  is  it  worthwhile  to  create  a  tropical  cyclone  model  which  uses 
potential  vorticity  as  a  predictive  variable?  One  study  which  suggests  the  answer  is  yes  is 
Haynes  and  McIntyre  (1987).  They  present  a  theorem  on  potential  vorticity  in  isentropic 
layers  which  they  state  as 

(i)  There  can  be  no  net  transport  of  Rossby-Ertel  potential  vorticity  (PV)  across  any 
isentropic  surface. 

(ii)  PV  can  neither  be  created  nor  destroyed,  within  a  layer  bounded  by  two  isentropic 
surfaces. 

This  holds  whether  or  not  there  is  diabatic  heating,  friction,  or  other  forces.  An  application 
of  this  theorem,  presented  as  a  thought  experiment,  leads  to  some  insight  into  tropical 
cyclones.  Their  experiment  involves  applying  a  pulse  of  cooling  to  an  unbounded,  stably 
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stratified,  rotating  fluid  at  rest  on  an  /-plane.  If  we  reverse  this  and  instead  apply  a  pulse 
of  heating  (Fig.  1.1)  we  have  a  system  analogous  to  a  tropical  cyclone,  the  heating  being 
the  release  of  latent  heat.  If  we  assume  the  heating  occurs  on  a  time  scale  much  faster 
than  any  dynamical  adjustment  process  then  by  the  arguments  of  Haynes  and  McIntyre 
the  isentropes  adjust  by  curving  as  shown  in  Fig.  1.2.  The  greatest  compaction  of  the 
isentropes  is  just  below  the  area  of  maximum  heating  and  this  creates  a  positive  potential 
vorticity  anomaly  in  that  region  while  the  reverse  is  true  above.  They  show  that  these 
anomalies  are  the  result  of  horizontal  transport  and  are  merely  a  result  of  redistribution 
of  potential  vorticity  along  the  isentropes.  Thus,  we  can  see  that  the  release  of  latent  heat 
leads  to  the  production  of  cyclonic  relative  vorticity  at  low  levels  and  anticyclonic  relative 
vorticity  at  upper  levels.  This  is  essentially  the  case  in  a  tropical  cyclone. 

Let  us  now  look  at  combining  the  potential  radius  with  the  potential  vorticity  to 
study  tropical  cyclones.  The  usefulness  of  this  technique  was  demonstrated  by  Schubert 
and  Alworth  (1987).  They  also  analyzed  the  balanced  vortex  equations  of  Eliassen  (1952). 
They  used  potential  radius,  as  the  horizontal  coordinate  with  potential  temperature  as 
the  vertical  coordinate.  This  combination  of  coordinates  led  to  a  simple  flux  form  of  the 
inverse  potential  vorticity  equation  which  separated  it  from  the  transverse  circulation. 
Solving  the  invertibility,  a  second-order  partial  differential  equation,  allowed  them  to  re¬ 
cover  the  wind  field  from  the  potential  vorticity.  The  solutions  show  that  latent  heat 
release  generates  potential  vorticity  at  low  levels  and  destroys  it  at  upper  levels.  In  addi¬ 
tion,  vertical  advection  causes  deepening  of  the  low  level  potential  vorticity  maximum  and 
pinching  off  of  the  upper  level  potential  vorticity  minimum.  Ultimately,  this  shows  that 
the  potential  vorticity  used  with  the  potential  radius  coordinate  makes  modeling  tropical 
cyclones  simpler  and  leads  to  better  physical  understanding. 

The  work  discussed  to  this  point  has  dealt  with  balanced  vortex  models  similar  to 
tropical  '•yclones,  but  not  a  closed  model.  Our  goal  is  to  move  toward  a  closed  model  which 
incorporates  these  principles.  We  begin  in  chapter  2  by  presenting  the  basic  concepts  of 
potential  vorticity  and  potential  radius  in  a  shallow  water  vortex.  We  start  with  the 
equations  of  motion  for  axisymmetric,  frictionless  flow  in  a  shallow  layer  of  homogeneous 


Figure  1.1:  Initial  state  of  experiment  proposed  by  Haynes  and  McIntyre  (1987).  A  pulse 
of  heating  (dashed  circle)  is  applied  to  an  unbounded,  stably  stratified,  rotating  fluid  at 
rest  on  an  /-plane. 
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Figure  1.2:  Final  state  of  experiment  proposed  by  Haynes  and  McIntyre  (1987).  Isentropes 
have  curved  in  response  to  heating  creating  positive  potential  vorticity  anomaly  (+)  below 
heating  and  negative  potential  vorticity  anomaly  (-)  above. 


incompressible  fluid  on  an  /-plane.  Then  we  introduce  the  potential  thickness,  the  inverse 
of  potential  vorticity,  and  write  the  equation  for  the  evolution  of  it.  This  equation,  in  some 
sense,  describes  the  release  of  latent  heat.  Next  we  transform  the  equations  to  potential 
radius  space.  This  eliminates  the  radial  wind  component  from  the  problem.  In  chapter  3 
we  present  several  different  ways  to  derive  and  solve  numerically  the  invertibility  principle. 
This  problem  turned  out  to  be  more  difficult  than  expected.  We  derived  several  different 
methods  solving  for  different  variables  with  three  different  numerical  methods-the  shooting 
method,  a  nonlinear  equation  solver,  and  a  tridiagonal  matrix  equation  solver.  Chapter  4 
gives  a  brief  description  of  Ooyama’s  model,  concentrating  on  the  prediction  of  the  wind 
field.  Then  in  chapter  5  we  generalize  these  ideas  to  the  full  tropical  cyclone  model.  The 
results  of  all  the  methods  are  presented  in  chapter  6. 


Chapter  2 


SHALLOW  WATER  MODEL 


We  start  by  presenting  the  concepts  of  potential  thickness,  potential  vorticity,  and 
the  invertability  principle  in  a  shallow  water  model. 

2.1  Potential  thickness  equation 

We  consider  axisymmetric,  frictionless  flow  in  a  shallow  layer  of  homogeneous  incom¬ 
pressible  fluid  on  an  /-plane.  If  we  assume  this  flow  is  in  gradient  balance,  the  gradient 
wind,  absolute  angular  momentum,  and  mass  continuity  equations  are 


('*;)•-& 

(2.1) 

(2.2) 

(2.3) 

d  d  d 

Dt  dt  +  Udr 

(2.4) 

where 


is  the  total  derivative,  u  and  v  are  the  radial  and  tangential  components  of  the  wind, 
<f>  =  g(h  -  h)  is  the  deviation  geopotential,  h  is  the  fluid  depth,  and  h  is  the  undisturbed 
depth  at  large  radius. 

The  vorticity  equation  for  this  system  is  derived  by  expanding  the  total  derivative  in 
(2.2)  to  obtain 

|w  +  «|;(n>)  +  ^|:(>a)  =  o. 

Now  multiplying  the  second  term  on  the  left  hand  side  by  r/r  and  operating  on  the  entire 
equation  with  d/(rdr)  yields 
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where  £  =  /  +  d(rv) / (rdr)  is  the  absolute  vorticity.  Now  we  use  this  equation  to  eliminate 
the  divergence  term  from  (2.3)  by  multiplying  (2.5)  by  -h  and  adding  C  times  (2.3)  to 


obtain 


^  Dt  h  Dt 


=  <Q- 


Defining  the  potential  thickness  H  =  (/ K)h  we  can  write 


DH 

Dt 


(2.6) 


Thus,  the  potential  thickness  is  the  depth  a  column  of  fluid  would  obtain  if  (  were  changed 
to  /  while  conserving  the  ratio  h/(.  It  is  important  to  note  the  potential  thickness  is 
inversely  proportional  to  the  potential  vorticity  £/h. 


2.2  Potential  radius  coordinate 


Let  us  define  the  potential  radius  R  by  fR2/2  =  ru  +  / r2/2,  the  right  side  of  which  is 
the  absolute  angular  momentum.  According  to  (2.2),  R  is  conserved  following  the  motion. 
We  now  define  T  =  t  and  transform  from  (r,  t )  space  to  {R,T)  space,  but  note  d/dt  implies 
fixed  r  while  d/dT  implies  fixed  R.  We  transform  the  partial  derivatives  so  that 


d_  _  dRd_  d_ 
dt  ~  dt  dR  +  dT' 

d_  _  d]ld_ 

dr  dr  dR  ’ 

The  definitions  of  R  and  (  allow  us  to  wrtite  (2.8)  as 


(2.7) 

(2.8) 


A  =  c  d 

rdr  f  RdR 


(2.9) 


Notice  in  regions  where  £  is  large  the  spatial  derivatives  are  stretched  in  R  space  giving 
better  resolution.  Next  we  transform  the  total  derivative  operator  (2.4)  by  substituting  in 
(2.7)  and  (2.8).  Then  using  the  definition  of  potential  radius  and  noting  the  conservation 
of  absolute  angular  momentum,  (2.4)  becomes 


D_  _  d_ 
Dt  ~  ar’ 


(2.10) 
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so  that  the  total  derivative  looks  like  a  local  derivative  in  ( R,T )  space.  This  allows  the 
potential  thickness  equation  (2.6)  to  be  written 


dH 

dT 


(2.11) 


Thus,  transforming  to  ( R,T )  space  allows  us  to  compute  the  potential  thickness  without 
explicitly  determining  the  radial  component  of  the  wind  u.  This  then  is  the  prognostic 
equation  of  the  model.  If  the  time  evolution  of  the  H  field  can  be  forecast  with  (2.11),  we 
can  use  it  to  solve  an  invertibility  principle.  The  derivation  and  solution  of  the  invertibility 
is  the  subject  of  chapter  three. 


Chapter  3 


THE  INVERTIBILITY  PRINCIPLE 


We  now  wish  to  derive  an  expression,  the  invertibility  principle,  which  enables  us  to 
recover  the  mass  field  h  form  the  potential  thickness  H.  The  invertibility  principle  for  this 
problem  can  be  stated  in  a  number  of  different  ways.  The  challenge  is  to  find  one  that  can 
be  solved  efficiently  on  a  personal  computer  and  easily  be  generalized  to  two  layers.  We 
have  explored  five  different  methods  of  solution  for  the  shallow  water  case  each  with  its 
own  strengths  and  weaknesses.  In  this  chapter  we  will  examine  these,  first  revealing  the 
various  forms  of  the  invertibility  and  then  discussing  the  numerical  integration  methods 
used  with  each. 

3.1  Method  one 


Method  one  involves  solving  for  a  new  variable  $. 
variables  $  and  V,  defined  by 

$  =  <£+  ;~v2, 

RV  =  tv. 


Let  us  introduce  the  two  new 

(3.1) 

(3.2) 


Substituting  (3.2)  into  the  definition  of  absolute  vorticity  and  transforming  the  result  into 
R  space  with  (2.9)  gives 


/  +  %? _ L 


d(RV)  ’ 


f  f~ 

while  from  (3.2)  and  the  definition  of  potential  radius  we  obtain 


(3.3) 


/+?  / 

t  ~  f  21'  ’ 

/  /  7T 


(3.4) 
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so  that  as  d(RV)/(RdR)  approaches  /  the  absolute  vorticity  becomes  infinite,  and  as  2 V/R 
approaches  /  the  absolute  circulation  per  unit  area  becomes  infinite.  The  definitions  (3.1), 
(3.2),  and  R  along  with  (3.3)  allow  us  to  transform  the  gradient  wind  equation  (2.1)  to 


(A 


We  now  wish  to  derive  the  relation  between  the  potential  function  $  and  the  potential 
thickness  H .  From  (3.1)  and  the  definition  of  potential  radius  we  get 

=  •  (3-6) 


Using  (2.9)  to  operate  on  r  and  using  the  definition  of  potential  thickness  we  get 

dr  _  R<f>* 
dR~  r<f>  ’ 


where  4>m  =  gH .  Using  the  definition  of  potential  radius  and  (3.2)  in  (3.5)  yields 

f2R(R 2  \ 


(£-')• 


We  solve  the  system  (3.6)-(3.8)  for  $  using  appropriate  boundary  conditions. 

At  the  center  of  the  vortex  symmetry  requires  the  inner  boundary  condition  to  be 

—  =  0  at  R  =  0.  (3.! 

UXt 


The  lateral  boundary  condition  $  — »  $  as  r  -+  oo  must  be  replaced  by  an  approximate 
condition  at  large  but  finite  r  =  R  =  Rf,.  To  derive  this  approximate  condition  suppose 
that  at  large  r  the  potential  thickness  has  not  been  disturbed  from  its  initial  value  so  that 
H  =  h  for  large  r.  We  also  assume  that  the  lateral  boundary  is  far  enough  away  from  any 
anomalies  in  the  potential  thickness  that  the  induced  rotational  flow  at  the  boundary  is 
weak  and  approximately  geostrophic.  Then,  the  weak  flow  conditions  |C  -  /I  <<  C  and 
-  $|  <<  $  allow  the  definition  of  potential  thickness  to  be  linearized  to  (//()-  1  + 
(4>  -  4>')/<t>‘  =  0,  (3.5)  to  fV  =  d$/dR,  and  (3.6)  to  $  =  <f>.  These  relations,  along  with 
(3.3),  can  be  combined  to  yield 


asO'S)-**-*’-0- 
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where  /u2  =  /2/$  and  $  is  the  undisturbed  or  basic  state  value.  The  solution  of  this  is 


$  =  AIo(tiRb)  +  BK0(nRb)  +  $, 


where  J0  and  Kq  are  the  modified  Bessel  functions.  We  eliminate  the  Jo  part  of  the 
solution,  because  $  must  be  bounded  as  R  — ►  oo  and  use  the  derivative  relation  for  A'o  to 


obtain 


where 


—  +  A($  ~  $)  =  0  at  large  J2, 
dK 


(3.10) 


(3.11) 


KofaRb)  '  v  ' 

Using  /  =  5.0  x  10~s  s-1,$  =  5.0  X  103  m2  s-1,  and  R b  =  1600  km  gives  A  =  9.77  x  10~7 
m-1.  Thus,  the  invertibility  principle  is  defined  by  (3.6)— (3.1 1 ). 

3.2  Methods  two  and  three 


Methods  two  and  three  both  solve  for  k  and  differ  only  in  the  numerical  schemes  used. 
We  use  the  definitions  of  potential  radius  and  the  deviation  geopotential  to  write  (2.1)  as 


4  \  r4  )  9rdr' 


Now  using  (2.9)  and  the  definition  of  potential  thickness  we  get 

dh  f2  HR  ( RA  A 


oh  __  f2  HR  f  R*  \ 
dR~  4g  h  V  r4  l) 


(3.12) 


which  relates  the  slope  of  the  top  boundary  to  the  wind  field,  through  the  difference 
between  the  potential  radius  and  the  actual  radius,  and  the  thermal  field,  through  the 
fluid  depth.  To  determine  r(R)  we  use  a  slightly  modified  version  of  (3.7) 


dr  RH 
dR  rh 


(3.13) 


The  boundary  conditions  here  are  the  same  as  in  method  one  with  h  replacing  $  in  (3.9) 
and  (3.10). 
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where  Rb  is  some  large  but  finite  R. 
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3.4  Method  five 


Method  five  proceeds  along  as  method  four  does,  but  continues  to  eliminate  h  and  r 
from  the  invertibility  equation.  By  using  the  relations  (3.3), 


ll  =(i-™)2 

R4  V  fRJ  ’ 


R2  +  r2  =  2  R2 


(-7.) 


(3.17a)  is  written 


d  (d(RV)\  1  dH  (  d{RV)\  V{f  *)(/  ®i) 

dR\RdR)  H  dRV  RdRj+  „(  t  2V\2 


0.  (3.18) 


which  becomes  the  new  statement  of  the  invertibility.  The  boundary  conditions  here  are 
the  same  as  in  method  four. 


3.5  Numerical  integration 


We  have  devised  three  different  numerical  analysis  schemes  to  solve  these  five  systems 
of  equations.  Methods  one  and  two  are  solved  using  the  shooting  method  with  the  4th 
order  Runge-Kutta  technique  to  determine  the  values  of  our  functions  at  each  grid  point 
and  the  secant  method  to  solve  the  lateral  boundary  condition,  methods  three  and  five 
using  a  nonlinear  equation  solver,  and  method  four  using  a  tridiagonal  matrix  equation 


solver. 


In  each  case  we  start  by  creating  a  grid  in  R  space  and  specifying  the  heating.  The 


grid  is  set  up  by  letting 


Rj  —  j  AfZ,  j  —  0, 1,2  ,...,</. 


(3.m 


In  addition  we  include  the  half  points 


#j+i  =  (j  +  ^AR,  j  =  0,1,2,...,  J 


(3.20) 
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where  A R  is  the  space  increment  in  the  horizontal.  The  heating  is  determined  from  (2,11). 
For  the  sake  of  illustrating  the  nature  of  the  solutions  of  this  equation  let  Q/h  =  S,  which 
we  assume  is  a  known  function  of  R  and  T.  The  solution  of  (2.11)  can  then  be  written 

H(R,T)  =  H(R,  0)exp  {jT 

An  interesting  consequence  of  (3.21)  is  that  it  is  not  necessary  to  compute  the  solution 
for  all  times  less  than  T  in  order  to  obtain  the  solution  at  time  T ;  only  the  “total  forcing” 
(the  exponent  in  (3.21))  need  be  specified.  If  two  solutions  happen  to  have  the  same 
distributions  of  potential  thickness,  they  also  have  the  same  distributions  of  mass  and 
wind  since  these  fields  are  derived  diagnostically  from  the  potential  thickness.  So  let  us 
consider  the  case  in  which  S  has  a  Gaussian  distribution  in  R  and  is  independent  of  time. 
Equation  (3.21)  then  reduces  to 

H{R,T)  =  H{R,  0)er<R\  (3.22) 


S(R,T')dr\. 


(3.21) 
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Figure  3.1:  The  potential  thickness  H  times  g  as  a  function  of  R  for  T  =  1,2, . . . ,  7  days. 
As  fluid  is  pumped  out  of  the  vortex,  the  potential  thickness  at  the  center  decreases. 
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3.5.1  The  shooting  method 

The  shooting  method  (Burden  and  Faires,  1985)  works  by  integrating  (or  shooting) 
from  the  center  to  the  outer  boundary  and  using  the  error  in  the  outer  boundary  condition 
to  make  another  pass.  The  values  of  the  functions  can  be  calculated  by  many  techniques. 
We  have  chosen  to  use  the  4th  order  Runge-Kutta  method.  We  use  the  secant  method  to 
reduce  the  error  in  the  outer  boundary  condition.  Our  algorithm  proceeds  as  follows.  We 
make  two  initial  guesses  for  the  variable  in  question.  Then,  stepping  outward,  calculate 
the  value  at  the  lateral  boundary  and  compute  an  error.  The  secant  method  is  then  used 
to  find  a  new  of  value  of  this  variable  which  is  used  to  start  another  iteration.  Convergence 
is  reached  when  the  error  is  reduced  to  an  acceptable  value.  The  4th  order  Runge-Kutta 
system  is  of  the  form 


(3.33) 
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or  more  specifically  in  the  case  of  method  three  F(/i(0))  =  0,  a  system  of  one  equation  in 
one  unknown,  and  in  method  five  F{Vj)  =  0,  a  system  of  79  equations  in  79  unknowns. 
Since  this  is  a  fairly  common  problem  we  decided  to  search  the  literature  for  existing 
software  to  use  as  a  “black  box”  to  solve  it.  We  found  a  suitable  subroutine  in  Kahaner 
et.  al.  (198S).  This  code  requires  us  to  provide  an  initial  guess  for  the  solution  and  a 
subroutine  to  calculate  the  value  '  the  function.  For  method  three  the  subroutine  that 
we  have  devised  uses  the  fourth-order  Runge-Kutta  method,  discussed  previously,  to  solve 
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for  h  and  r  with  an  initial  guess  of  h( 0)  =  ha-  For  method  five  the  subroutine  solves  the 
finite  difference  form  of  (3.18) 

1  (Rj+iVj+i  -  RjVj  RjVj  -  \ 

J  (Ai?)2  v  rj+l  rh  r 

L- ««•> 

with  initial  guesses  of  V,  =  0. 

3.5.3  The  tridiagonal  matrix  solver 

When  we  write  (3.17a)  in  finite  difference  form  we  create  a  system  of  linear  equations 
in  V .  If  we  use  matrix  notation  to  solve  this  system,  the  coefficient  matrix  of  V  is 
tridiagonal,  all  the  entries  are  zero  except  those  on  the  diagonal  and  on  both  sides  of  it. 
This  form  is  particularly  convenient  when  solving  the  matrix  equation 


Ax  =  b. 


(3.40) 


One  method  to  accomplish  this  is  Crout  reduction(Burden  and  Faires,  1985).  If  A  is  the 
coefficient  matrix 


A  = 


/  Oil  ®12  0 

021  022  <*23 

0  **. 


0  \ 


an-l,n 


\  0  ...  0  On,n- 1  ®nn  J 

it  can  be  factored  into  the  lower-triangular  matrix  L  and  the  upper-triangular  matrix  U 
(LU-decomposition) 


L  = 


(hi  0 
hi  I22 

0  **. 

•  • 

•  • 

V  0  ... 


0  \ 


0  1 


n,n- 1  ^nn  ' 
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Since  A  =  LU,  (3.40)  becomes  LUx  =  b  or 

Ux  =  L-1b 

which  we  use  to  solve  for  x.  The  LU-decomposition  is  carried  out  by  noting  the  multipli¬ 
cations  involved  in  A  =  LU 

an  =  hi, 

ais-i  =  /M_ i  for  each  t  =  2, 3, . . . ,  n, 

a„  =  /«,,•_ i u,_ 1 ,i  +  hi  for  each  t  =  2,3,...,n, 

=  /„•«,,,+ 1  for  each  *  =  1,2,. . .  ,n  -  1. 

Once  the  matrices  L  and  U  are  known  we  introduce  a  dummy  variable  z  and  solve 

Lz  =  b 

using 


)  for  each  i  =  2, 3, . . . ,  n. 

The  final  step  is  to  solve 

Ux  =  z 

using 

xn  —  zn 

=  2j  —  u,il+iZj+i  for  each  t  =  n-l,...,l. 

This  form  is  particularly  attractive  in  our  problem,  because  it  updates  all  of  the  values  at 
the  interior  grid  points  at  one  time  without  changing  the  values  at  the  end  points  which 
in  our  case  are  always  zero. 
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So  we  write  (3.17a)  in  finite  difference  form 

gh3  (  R:+1VJ+1  -  R:V3  RM-Ri- iVW 


/2(A£)2 


f  Vj+1  -  RjVj  RjVj  -  Rj-xVj-x  \  r  R](R]  ±  rj)  H]  \ 
1  Rj+i  -  R}.k  J  +  l  2r1  h2J 

—  a( 

f\dRjj 


(3.41) 


or  the  more  convenient  form 


(  R]~l  \  y  1 

f  2  R]  R](R*,+rl)n}fH&R)*) 

[y.  _  (Rj+y) 

U.-.J  '-,  +  i 

j  2 

IR,-lR,+i  2r<  A?  gh,  j 

<Vl  Uj+J 

Next  we  determine  the  coefficient  matrix  A  and  the  right  hand  side  of  (3.42).  These 
are  passed  to  a  subroutine  that  performs  the  LU-decomposition  and  solves  the  tridiagonal 
system.  The  subroutine  passes  back  a  solution  for  V  which  we  use  to  update  r2  and  h 
and  compute  a  residual. 


Chapter  4 


REVIEW  OF  THE  OOYAMA  TROPICAL  CYCLONE  MODEL 

4.1  The  fluid  system  and  equations  of  motion 

The  cyclone-scale  motions  in  Ooyama’s  model  are  approximated  by  a  quasi-balanced 
axisymmetric  vortex.  It  has  three  vertical  levels  two  of  which  are  prognostic;  the  other 
is  diagnostic.  For  simplicity  the  fluid  is  assumed  to  be  homogeneous  and  incompressible. 
This  assumption  eliminates  any  true  thermodynamics  which  is  reasonable  for  adiabatic 
motions,  but  can  be  a  problem  for  diabatic  processes.  Fortunately,  only  unstable  moist 
convection  involves  diabatic  motions  and  this  is  formulated  indirectly  in  the  model. 

The  system  consists  of  two  main  layers  and  a  shallow  planetary  boundary  layer.  Fig. 
4.1  shows  a  schematic  diagram  of  the  fluid  system.  This  allows  for  the  low  level  cyclone 
and  the  upper  level  anticyclone  of  the  tropical  storm.  The  density  of  each  layer  is  assumed 
constant,  and  the  system  is  stably  stratified  such  that 

e  =  —  <  1  and  p0  =  px.  (4.1) 

P  i 

Also  important  to  note  are  the  boundary  layer  has  constant  thickness  ho,  and  the  top 
boundary  of  the  system  is  a  free  surface. 

For  convenience  he  introduces  a  reference  density  p  such  that 

P  =  P\  =  Po  and  p2  =  ep. 

He  also  defines  the  nondimensional  static  stability 

t r=l-£ .  (4.2) 

Since  the  fluid  is  assumed  to  be  incompressible  changes  in  the  thickness  correspond 
to  adiabatic  temperature  changes.  A  layer’s  thickness  can  only  change  by  convergence 


Figure  4.1:  Schematic  diagram  of  the  fluid  system  proposed  by  Ooyama  (1969). 
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and  divergence  of  the  fluid  caused  by  the  radial  wind.  The  tangential  wind  can  not  have 
any  effect  on  the  divergence,  because  of  axisymmetry.  Thus  the  radial  mass  fluxes,  for 


each  layer  are 


tpo  =  -honor, 


V»i  =  -hyuir, 


■02  =  -eh2u  ir,  (4.5) 

where  Uj  is  the  radial  component  of  the  wind  defined  positive  outward.  Actually  the  mass 
flux  is  pty,  but  tfjj  is  used  with  little  confusion. 

Diabatic  processes  are  represented  by  the  change  in  density  of  a  fluid  parcel  associated 
with  the  vertical  motion.  Thus,  diabatic  heating,  the  movement  of  air  from  high  density 
to  lower  density,  is  pQ+ ,  and  diabatic  cooling  is  pQ~  where  Q+  and  Q~  are  the  diabatic 
mass  fluxes. 

Next  he  writes  the  mass  continuity  equations.  Since  the  boundary  layer  has  constant 


thickness  its  continuity  equation  is 


n  d^° 

0  =  —  -  tn 

TOT 


where  w  is  the  vertical  velocity  at  the  top  of  the  layer.  The  equations  for  the  main  layers 


dhi  dxpi 

JT^7^-Q  +  W 


dh2  0t/>2  _ 

c-dT  =  7s;+Q 


where  Q  =  Q+  -  Q  is  the  net  diabatic  flux. 

Using  the  hydrostatic  approximation  the  pressure  in  each  layer  is  written 


Pi  =  <7p(*o  +  hi  +  eh2  —  z)i 


Pi  =  £SP(^o  +  hi  +  h2  -  z), 


(4.10) 
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Po  =  Pi 


(4.11) 


where  z  is  the  height  above  the  sea  level  and  g  is  the  acceleration  due  to  gravity.  Now, 
if  hj  and  pj  denote  the  undisturbed  or  average  values  of  hj  and  pj  then  the  deviation 


geopotential  is  defined  by 


*  = 

Pj 


(4.12) 


As  previously  stated  the  tangential  circulation  is  represented  by  a  balanced  axisym- 
metric  vortex.  The  radial  equation  of  motion  is  thus  approximated  by  the  gradient  wind 


equation 


(/+?),= 


(4.13) 


where  /  is  assumed  constant.  Since  from  (4.12)  <f> o  —  4>\,  from  (4.13)  uo  =  t>i- 
Ooyama  defines  the  absolute  angular  momentum  for  this  system  as 

Mj  =  VjT  +  i/r2.  (4.14) 

The  angular  momentum  budgets  in  each  of  the  main  layers  is  then  written 

Q  A 

^(AiA/j)  =  ^(0iM,  +  Ai)  -  Zu  +  Z0,  (4.15) 


0  d 

£  =  +  A2)  +  Z12 


(4.16) 


where  A j  are  the  radial  fluxes  of  angular  momentum  due  to  lateral  eddy  transport,  Z 12 
is  the  vertical  flux  of  angular  momentum  between  the  two  main  layers,  and  Zo\  is  the 
vertical  flux  of  angular  momentum  from  the  boundary  layer  to  the  lower  main  layer.  The 
model  equations  are  derived  from  these  basic  equations. 

4.2  The  model  equations 


The  equations  for  Ooyama’s  model  are  listed  in  the  box  below.  A  model  run  consists  of 
predicting  <f>\  and  $2  with  (4.17)  and  (4.18).  In  these  equations  Gi  are  the  diabatic  heating 
terms  which  are  dependent  on  Q+.  This  in  turn  depends  on  77  which  is  a  nondimensional 
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parameter  developed  from  energy  considerations  and  is  dependent  on  representative  values 
in  layers  0  and  1  of  9e ,  an  arbitrary  constant  0,  and  the  saturation  value  in  layer  2  0*. 
Once  4>i  is  known  hi  and  /12  are  determined  with  (4.19)  and  (4.20)  after  which  V\  and  V2 
are  predicted  with  (4.21)  and  (4.22)  where  F/  are  the  friction  terms. 


Chapter  5 


INCLUSION  OF  POTENTIAL  VORTICITY  IN  THE  OOYAMA 
TROPICAL  CYCLONE  MODEL 

The  dosed  tropical  cyclone  model  is  developed  from  the  basic  concepts  of  the  shallow 
water  model  with  some  major  modifications.  We  use  a  more  general  fluid  system  with 
three  layers.  The  lowest  layer  represents  the  moist,  cyclonic  boundary  layer,  and  the  upper 
two  layers  represent  the  cyclonic  and  anticyclonic  flows  above.  In  addition  to  the  more 
complicated  system  we  include  surface  and  internal  frictional  forces.  Absolute  angular 
momentum  is  no  longer  conserved  with  friction  included.  However,  using  the  potential 
radius  coordinate  remains  an  advantage.  In  this  chapter  we  generalize  the  concepts  of 
the  shallow  water  model  and  derive  and  solve  the  invertibility  principle  for  the  two  main 
layers  of  the  tropical  cyclone  model. 

5.1  Potential  thickness  equations 

We  use  the  fluid  system  proposed  by  Ooyama  consisting  of  three  layers,  the  lowest 
two  having  density  p,  and  the  uppermost  having  density  ep  where  e  <  1  so  that  the  system 
is  statically  stable.  The  layer  thicknesses  are  ho,  h\,  and  h?  where  ho  is  a  constant,  and 
hy  and  /12  depend  on  radius  and  time.  We  define  the  deviation  geopotentials 

<t>\  =  s[(hi  -  hi )  +  e(h2  -  h2)]  (5.1 ) 

and 

4>2  =  p[(hi  -  hi)  +  (/12  -  /12)]  (5.2) 

and  write  the  gradient  wind  equations  for  the  upper  two  layers  as 


29 


where  l  is  the  layer  number.  Since  we  assume  the  hydrostatic  approximation  and  the  lower 
two  layers  have  the  same  density,  they  also  have  the  same  radial  pressure  gradient  force. 
If  we  further  assume  that  the  boundary  layer  is  also  in  gradient  wind  balance,  we  conclude 
that  vq  =  v\.  Although  there  is  no  difference  between  the  tangential  wind  components  in 
layers  0  and  1,  we  shall  allow  a  difference  between  their  radial  components. 

We  now  turn  our  attention  to  the  mass  conservation  and  angular  momentum  princi¬ 
ples,  both  of  which  provide  prognostic  equations  for  the  layers  above  the  boundary  layer. 
An  objection  sometimes  raised  against  the  use  of  a  homogeneous  incompressible  fluid  sys¬ 
tem  for  tropical  cyclone  modeling  is  that  it  does  not  include  explicit  thermodynamics. 
This  objection  can  be  answered  by  means  of  analogy.  Suppose  we  consider  constant  po¬ 
tential  temperature  layers  instead  of  constant  density  layers.  Then,  the  geopotential  is 
replaced  by  the  Montgomery  potential  and  the  layer  thicknesses  are  replaced  by  the  pres¬ 
sure  differences  across  the  isentropic  layers.  Except  for  some  non-Boussinesq  effects  the 
isentropic  layer  equations  turn  out  to  be  identical  to  the  isopycnic  layer  equations.  With 
this  in  mind,  let  us  consider  the  isopycnic  layer  continuity  equations 

Q 

( ru0h0 )  +  w  =  0,  (5.4) 

1ft  +  rfr  (rUlhl)  +  ~  ~  w  =  °>  (5-5) 

£1F  + jJr(erU2/l2  )~9+  +  <2~  =  0’  (5-6) 

where  w  is  the  vertical  velocity  at  the  top  of  the  boundary  layer,  pQ+  is  a  vertical  mass 
flux  from  layer  1  to  layer  2,  and  pQ~  is  a  vertical  mass  flux  from  layer  2  to  layer  1.  In 
the  isentropic  layer  model  the  analogue  of  the  net  mass  flux  p(Q+  —  Q~ )  is  06,  where 
a  =  —dp/d6  is  the  isentropic  pseudo-density.  Thus,  Q+  represents  movement  of  air  from 
low  to  high  potential  temperature  (diabatic  heating)  while  Q~  represents  movement  in 
the  opposite  sense  (diabatic  cooling). 

Let  us  define  the  absolute  angular  momentum  m\  and  the  potential  radius  Ri  as 
mi  =  fltf/2  =  rvi  +  /r2/ 2.  Then,  the  absolute  angular  momentum  equations  can  be 
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written 

'  £cD  |vi|  t>i,  /  =  0 

fRiRi=  <  j^{Q~  +m)(«2  ~  «i)>  J=1  (5.7) 

.  +  p)(v 1  -  ”2),  *  =  2 

where  c#  is  the  drag  coefficient,  /i  is  the  constant  coefficient  for  momentum  transfer 
through  shearing  stress  at  the  interface,  and  Ri  is  the  total  derivative  of  Ri.  We  have 
neglected  lateral  momentum  transport  due  to  turbulent  eddy  processes. 

Assuming  Q+  and  Q~  can  somehow  be  parameterized,  we  regard  (5.1)-{5.7)  as  a 
closed  system  in  uo,  tti,  u2,  t>i,  u2,  tu,  fa,  fa,  hi,  h2.  Since  (5.1)-(5.3)  diagnostically  relate 
h\,  fa,  Vi,  v2,  the  prediction  of  hi  and  hi  from  (5.5)  and  (5.6)  must  be  consistent  with 
the  prediction  of  vj  and  u2  from  (5.7).  This  implies  a  coupled  pair  of  diagnostic  equations 
for  uj  and  u2  (see  Equation  (4.1)  of  Ooyama,  1969).  This  coupled  pair  can  then  replace 
either  (5.5)  and  (5.6)  or  (5.7).  In  the  next  section  we  adopt  a  different  approach.  We  first 
combine  (5.5)  and  (5.6)  with  (5.7)  to  obtain  the  potential  thickness  (or  inverse  potential 
vorticity)  equations.  These  equations  contain  ut  and  u2  as  advecting  velocities.  We  then 
transform  the  independent  variable  from  actual  radius  to  potential  radius,  which  removes 
«i  and  u2  from  the  problem. 

Defining  the  vertical  component  of  absolute  vorticity  by 


we  can  differentiate  (5.7)  to  obtain 

§  +  4;  K'  -  /**)  =  <5-9> 

where  /RiRi  has  been  used  as  shorthand  notation  for  the  right  hand  sides  of  (5.7).  Defining 
the  potential  vorticity  as  qi  =  Q/hi  we  can  write  (5.9)  as 

Q  Q 

^(Mi)  +  ^  (ru ih,qi  -  fRiR^j  =  0.  (5.10) 

Thus,  ruihiqi  -  fRiRi  can  be  interpreted  as  the  flux  of  potential  vorticity.  This  flux  is 
along  the  layer  so  that  the  interfaces  between  layers  are  impermeable  to  potential  vorticity. 
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Integrating  over  the  area  within  some  large  radius  at  which  both  rui  and  RiR{  vanish  we 


obtain 


d  f 

« j  ,>}h‘Tdr  =  °- 


(5.11) 


In  this  sense  the  potential  vorticity  within  the  layer  is  indestructible.  These  two  restrictions 
on  the  way  in  which  the  potential  vorticity  can  change  have  been  clearly  discussed  by 
Haynes  and  McIntyre  (1987)  and  McIntyre  (1987). 

We  can  now  combine  (5.5)  and  (5.6)  with  (5.9)  to  obtain  the  potential  thickness 
equations.  We  do  this  by  multiplying  (5.5)  and  (5.6)  by  f/Cb  Now  we  subtract  (5.9) 
times  ( /hi)/Ci  fr°m  (5.5)  and  (5.9)  times  (tfhi)l(%  from  (5.6)  to  obtain 

i 


+  -^-(2+-Q--tn)  =  ° 


(5.12) 


where 


sh2  wr,  /  e  , 
~ST  +  "2~db  + 

.^(Q+_Q-)  =  0 


H,  =  yh, 

Q 


(5.13) 


(5.14) 


is  the  potential  thickness  of  layer  l.  Since  Hi  is  the  inverse  of  the  potential  vorticity,  we  can 
interpret  Hi  as  the  depth  layer  l  would  acquire  if  (/  were  changed  to  /  under  conservation 
of  the  ratio  h(/Q. 

5.2  Potential  radius  coordinate 


With  (r,  t)  as  independent  variables  we  regard  the  Ri(r,t)  as  dependent  variables. 
We  now  reverse  this  situation  and  regard  ( R,T )  as  independent  variables  and  ri(R,T) 
as  dependent  variables.  We  define  T  =  t  but  note  that  d/dt  implies  fixed  r  while  d/dT 
implies  fixed  R.  We  transform  the  partial  derivatives  to 

d  dRi  d  d_ 
dt  dt  dR  +  dT 


(5.15) 
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and 


d  dRi  d 
dr  dr  dR 


(5.16) 


the  second  of  which  combined  with  the  definitions  of  potential  radius  and  absolute  vorticity 


is  written 

d  _  o  a 

rdr  f  RdR' 

Using  (5.15)  and  (5.16)  we  write  the  total  derivative 


(5.17) 


a  ^  d  _  d  A 

dt  +  U,dr  dT  + 


(dRi  ,  dR,\  a 

[dt  +U‘  dr  )  dR 


a  ^  d  _  d  d 

dt+Uldr  dT  +  RldR ' 
Substituting  (5.17)  and  (5.1 8)  into  (5.12)  and  (5.13)  gives 


(5.18) 


(5.19) 


and 

(5-20) 

which  now  serve  as  the  prognostic  equations  of  the  model.  Note  that  the  radial  components 
u\  no  longer  appear  explicitly. 


5.3  Invertibility  principle 


Again  we  wish  to  derive  expressions  for  the  relationships  between  hi  and  Hi.  Following 
the  procedures  of  method  five  for  the  shallow  water  model  we  start  with  the  definition  of 
potential  thickness  (5.14)  and  then  introduce  the  new  variables  V)  defined  by 

RV,  =  rv,.  (5.21) 


Substituting  (5.21)  into  (5.8)  and  transforming  the  result  into  R  space  with  (5.17)  gives 


(5.22) 


while  from  (5.21)  and  the  definition  of  potential  radius  we  obtain 


f+27L  _  / 

/  -f-T 


(5.23) 
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so  that  as  the  absolute  vorticity  becomes  infinite  d(RV[)/(RdR)  approaches  /,  and  as  the 
absolute  circulation  per  unit  area  becomes  infinite  2 Vi/R  approaches  /  . 

Now  using  (5.22)  in  (5.14)  and  taking  d/dR  of  the  new  equation  gives 


dR[ 

<d(RV i)\  ffdht 

K  RdR  )  hi  Ci  dR 

(t  d(RV i)> 
V  RdR  t 

\  1  dH  i 

(5.24) 

and 

9  1 
dR' 

<d(RV2)\  f  fdh2 

V  RdR  )  '  h2( 2  dR 

(f  d{RV2 )> 
V  RdR  j 

v  1 

'  H2  dR 

(5.25) 

which  can  be  written 

_±l 

HR' 

'd(RV i)\ 
v  RdR  ) 

1  dHi(  d(RV,)\ 
Hi  dR  V  RdR  ) 

1  +  crgHi  (o> 

\2  dhl  n 

)  "a  =0 

(5.26) 

and 

OR' 

< d(RV2)\ 
K  RdR  ) 

1  dH2f  d{RV2 )} 
H2  dR  V  RdR  ) 

1  +  crgH2{  (2/ 

\ 2  dh2 
)°SJR=0 

(5.27) 

where  a  s=  1  -  e.  Now  we  need  to  find  an  expression  for  dhi/dR.  Equations  (5.1),  (5.2), 
and  the  definition  of  potential  radius  are  used  to  write  (5.3)  as 

4  \  rj  /  rdr  ’ 

From  these  two  equations,  (5.1),  and  (5.2)  we  write 

f2(Rj-eR4  \  df± 

4  V  <7r4  )  9  rdr 
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(5.28).  To  calculate  hi  we  take  some  arbitrary  R.  Then  we  need  to  find  the  coresponding 


R2  to  subtract  from  it.  Because  of  the  coordinate  transformation  r  space  is  distorted  (Fig. 
5.1),  and  we  can’t  simply  use  the  value  of  R  in  layer  two.  We  have  to  find  the  value  of  rj 
coresponding  to  R  and  then  find  R2  on  this  r  surface.  This  is  depicted  graphically  in  Fig. 
5.2.  To  do  this  we  need  to  know  rj  as  a  function  of  R.  This  will  be  discussed  further  in 
the  section  on  the  numerical  integration. 

So,  we  substitute  (5.28)  and  (5.29)  into  (5.26)  and  (5.27)  to  get 

d  fd(RVx)\  _  1  dHx(  TO, 
dR\  RdR  )  Hi  dR  V  RdR  J 


M*)I 


(5.30a) 


d  fd(RV2)\  1  8H2 


IN  _  1  9H2 
)  H2  dR 
3 

R(R*-[1 i 
4  \  [r2 


d(RV2)' 


( f  afHValV 

V  RdR  )  RfR'-lRijr^R))}* 


ffl)  = 


(5.30  b) 


o9H2  4  \  [r2(R)Y 

which  along  with  appropriate  boundary  conditions  are  the  invertibility  for  the  two  layer 


model. 


At  the  center  of  the  vortex  symmetry  requires  that 


Vi  =  0  at  R  =  0, 


(5.30c) 


and  we  assume  the  lateral  boundary  is  far  enough  from  the  center  that  the  flow  decays  to 
zero.  So 

V,  =  0  at  R  =  Rb  (5.30d) 

where  R (,  is  some  large  but  finite  R. 

A  model  computational  cycle  consists  of  predicting  new  values  of  Hi  and  H2  using 
(5.19)  and  (5.20)  followed  by  solution  of  (5.30)  to  obtain  Vj.  We  then  solve  (5.21)  for  v /. 


R 


Figure  5.1:  The  transformed  coordinate  system.  Surfaces  of  r  =  constant  are  distorted  in 
R  space. 
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5.4  Numerical  integration 


We  generalize  the  nonlinear  equation  solver  method  developed  for  the  shallow  water 
problem  to  the  tropical  cyclone  model.  Our  problem  now  is  solving  the  system 

Wk)  =  o 

We  again  use  the  subroutine  from  Kahaner  et  al.(1988)  providing  guesses  Vji;  =  0.  Our 
system  now  is  158  equations  in  158  unknowns. 

The  subroutine  we  provide  solves  the  finite  difference  forms  of  (5.30a)  and  (5.30b). 
However,  to  do  so  we  need  to  know  r/(i2)  which  we  find  from  (5.21) 


-=4-7^ 

The  finite  difference  forms  of  (5.30a)  and  (5.30b)  are 

p  _  i  ( Rj+\V\,j+i  -  -RjVi.j  -  , 

,J"(A*n  r 

(±M±\  (f  Ri+xVij+x-Ri-xVij-i\ 

\Hl  dR  JjV  2 RjAR  ) 

V_ _ )  R;  ( R)  -  e[R2{rx{Rj))]A  \ 

vgHij  4  V  [rtiRj)}4  °) 


(5.31) 

(5.32) 


(5.33) 


and 


1  ( Rj+}V2j+i  -  RjV2j  RjV2j-Rj- iV2lj_i\ 

2,J‘(A  r)A  rh  r 

Rj+iVij+i  -  Rj-iXhj-i 


(  1  9H2\  (  f  Rj+iV2,j+\  -  Rj- lVij-A 
\H2dRj:V  2RjAR  ) 


( / 


L)  Rj  (R]  -  IRxMmV 


)- 


vglh.,  4  V  [r2(*,)]<  )  ”  (53J) 

The  algorithm  proceeds  as  in  the  single  layer  case  with  the  grid  set  up  in  R  space. 
We  then  calculate  r/j  at  the  grid  points  and  use  a  cubic  spline  routine  as  another  “black 
box”  to  create  a  smooth  function.  We  also  use  the  cubic  spline  to  create  R2(r\(R:))  and 
#i(r2(#j))-  Then  when  we  need  values  of  these  we  evaluate  the  splines. 
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To  calculate  the  heating  we  assume  frictionless  flow  with  w  =  0  and  Q  =  0.  This 


reduces  (5.19) to 


and  (5.20)  to 


fr=° 


®Jh.-Q+Ei  =  o. 

3T  V  eh2 


(5.35) 


(5.36) 


We  now  let  5/  =  /hi  be  some  known  function  of  R  and  T.  Now  we  integrate  (5.35) 
and  (5.36)  to  obtain 


H\(R,T)  =  Hi(Rt0)  exp  (-  j*  S1(R,T')dT'Sj 


(5.37) 


H2(R,T)  =  H2(R,0)  exp  S2(R,T')dT' 


(5.38) 


Now  to  illustrate  the  solutions  of  (5.37)  and  (5.38)  we  let  S  =  S\  =  S2/e  be  a  function 
with  a  Gaussian  distribution  in  R  and  independent  of  time.  Then  (5.37)  and  (5.38)  reduce 


H\(R,T)  =  Hi(R,  0)eTi{R) 


(5.39) 


H2(R,T )  =  H2(R,0)e^W 


(5.40) 


where  t\(R)  =  -SoTexp(-R2 /R%)  and  t2(R)  =  SqT exp(-R2  /  Rq).  This  creates  a  mass 
sink  in  layer  one  and  a  mass  source  in  layer  two. 

After  Hi  is  computed  at  all  the  grid  points,  the  solution  is  carried  out  as  in  the  single 
layer  case.  The  output  variables  U(,h(,  and  Q/f  are  calculated  from 


vi.i  = 


(5.41) 


h,J  ~  H,J  f  _  g(g.vj..)  ’ 


(5.42) 


Oj  h‘.j 

f  Hi/ 


(5.43) 


Chapter  6 


RESULTS 


6.1  Shallow  water  case 

Four  of  the  five  methods  used  to  solve  the  shallow  water  invertibility  work  very  well. 
They  are  efficient  and  able  to  produce  fairly  strong  vortices.  The  fifth  is  efficient,  but  can 
produce  only  a  weak  vortex.  For  comparison,  each  method  was  run  at  168  hours  except 
method  four.  Method  four  does  not  run  that  far,  so  it  was  run  at  96  hours.  Fig.  6.1  shows 
the  v  and  (//  fields  for  all  five  methods.  All  the  methods  produce  very  similar  results 
that  resemble  a  tropical  cyclone.  The  individual  results  are  now  presented  in  more  detail. 

6.1.1  Methods  one  and  two 

Method  one  works  very  well.  The  v  and  £//  fields  for  this  method  are  shown  in  Fig. 
6.2.  It  runs  to  a  maximum  of  242  hours.  At  168  hours  it  produces  a  peak  wind  of  36.66 
m/s  at  R  =  180  km.  The  error  in  the  boundary  condition  is  reduced  seven  orders  of 
magnitude  to  10"9  in  five  iterations. 

Method  two  also  works  very  well.  It  runs  to  a  maximum  of  330  hours  and  at  168 
hours  (Fig.  6.3)  produces  a  peak  wind  of  36.41  m/s  at  R  =  180  km.  The  error  in  the 
boundary  condition  here  is  reduced  seven  orders  of  magnitude  to  10-11  in  seven  iterations. 

However,  these  methods  are  hard  to  generalize  to  two  layers,  because  of  the  secant 
method.  The  secant  method  is  used  to  solve  a  function  of  a  single  variable.  The  two  layer 
model  when  formulated  as  in  methods  one  and  two  is  a  problem  of  two  equations  with 
two  unknowns.  Specifically,  the  errors  in  the  lateral  boundary  conditions  of  each  layer  are 
dependent  on  the  depths  of  both  layers  at  the  center.  Thus,  the  secant  method  can  not 
be  used  to  solve  this  problem. 


Figure  6.1:  Top:  The  tangential  wind  (u)  as  a  function  of  R  for  all  five  solution  methods 
and  T  =  7  days  (Method  four  T  =  4  days.).  Bottom:  The  normalized  vorticity  (C//)  as 
a  function  of  R.  All  else  the  same. 
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6.1.2  Method  three 

Method  three  is  slightly  slower  than  one  and  two,  but  still  works  very  well.  It  runs  to 
a  maximum  of  360  hours.  At  168  hours  (Fig.  6.4)  it  produces  a  peak  wind  of  36.41  m/s 
at  R  =  180  km.  The  error  in  the  lateral  boundary  condition  is  reduced  thirteen  orders 
of  magnitude  to  10-11  in  nine  iterations.  Again  this  method  is  hard  to  generalize.  The 
problem  is  not  well  defined  in  two  layers.  The  error  in  the  outer  boundary  condition  is  not 
uniquely  determined  by  the  guesses  of  h/(0).  The  problem  is  the  coordinate  transforma¬ 
tion.  As  previously  stated  when  we  make  the  coordinate  transformation  and  discretetize 
the  two  layer  model  we  have  to  interpolate  R/(rn(R)).  When  we  use  the  nonlinear  equation 
solver  the  interpolation  is  done  separately  and  is  not  constrained  in  any  way.  Thus,  if  we 
choose  different  interpolations  we  would  expect  the  nonlinear  equation  solver  to  produce 
different  results.  So,  we  can  not  uniquely  define  the  problem  to  be  solved. 

6.1.3  Method  four 

Method  four  does  not  work  well.  It  runs  to  a  maximum  of  only  108  hours  and  to  get 
it  that  far  we  have  to  use  previous  solutions  as  initial  guesses.  At  96  hours  (Fig.  6.5)  it 
produces  a  peak  wind  of  16.43  m/s  at  R  =  180  km.  The  problem  here  seems  to  be  in 
computing  the  diagonal  of  the  coefficient  matrix  for  (3.42).  The  tridiagonal  solver  gets 
very  sensitive  as  the  vorticity  gets  large  and  produces  infinite  vorticity  in  a  fairly  short 
time. 

6.1.4  Method  five 

Method  five  is  the  one  we’ve  been  looking  for.  It  works  well,  but  is  somewhat  slow. 
While  the  other  methods  all  ran  in  seconds,  this  one  runs  in  approximately  two  minutes. 
However,  we  were  able  to  generalize  it  with  good  results.  In  the  shallow  water  case  at  168 
hours  (Fig.  6.6)  it  produces  a  peak  wind  of  36.93  m/s  at  R  =  160  km. 

Since  this  method  is  the  one  we  used  to  solve  the  two  layer  case  a  more  complete 
output  set  is  included.  Fig.  6.7  shows  the  V(R)  fields  corresponding  to  the  H(R,T)  fields 
in  Fig.  3.1.  Fig.  6.8  shows  the  corresponding  t>(r)  fields  computed  from  (3.2).  Note  that 
as  the  r  at  which  the  peak  u(r)  occurs  decreases  with  time,  the  R  at  which  the  peak  V(R) 


zeta/f  v(m/s) 
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Figure  6.5:  Top:  The  tangential  wind  (t>)  as  a  function  of  R  for  method  four  (T  =  4  day 
Bottom:  The  normalized  vorticity  ((//)  as  a  function  of  R  for  method  four  (T  =  4  day 


0 


200 


400 


600 


600 


R(km) 

Figure  6.6:  Top:  The  tangential  wind  (v)  as  a  function  of  R  for  method  five  (T  =  7  days). 
Bottom:  The  normalized  vorticity  ((//)  as  a  function  of  R  for  method  five  (T  =  7  days). 
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occurs  increases  with  time.  This  is  a  consequence  of  (3.3) — stretching  in  regions  of  large 
vorticity. 

Another  way  to  view  the  solutions  is  shown  in  Figs.  6.9-6.11.  Fig.  6.9  shows  the  r 
space  distributions  of  H  and  Figs.  6.10  and  6.11  show  its  factors  h  and  / /(.  We  see  that 
as  the  mass  sink  forces  changes  in  the  potential  thickness,  the  fluid  responds  by  mutually 
adjusting  the  mass  field,  h,  and  the  wind  field,  (//,  to  maintain  gradient  balance.  Which 
is  adjusted  more,  the  mass  field  or  the  wind  field?  From  Fig.  6.9  H  at  the  center  of  the 
vortex  is  initially  fairly  low.  At  the  same  time  h  is  very  high  while  £//  is  small.  Thus, 
initially  the  wind  field  makes  a  greater  contribution  to  H.  As  the  system  evolves,  H  gets 
lower  as  does  h,  but  (/f  gets  very  large.  So,  with  time  the  contributions  from  h  and  (J  f 
approach  each  other.  Thus,  the  mass  field  adjusts  to  the  wind  field. 

6.2  The  two  layer  model 

The  results  shown  for  the  two  layer  model  are  for  the  heating  given  in  (5.39)  and 
(5.40)  with  tfi(.R,0)  =  5133.0  m,  H2(R, 0)  =  3523.8  m,  Ro  =  250  km,  and  r2  =  -0.5rj. 
The  model  was  run  at  SqT  =  1,2,..., 7  and  /  =  5.0  x  10_5s_1.  Fig.  6.12  shows  the 
H(R,T )  field  in  layer  one  for  the  times  listed  above. 

The  results  show  good  tropical  cyclone  structure.  The  lower  layer  is  as  in  the  single 
layer  case.  Fig.  6.13  shows  V(R)  and  again  we  get  the  radius  of  maximum  “wind” 
increasing  with  time  while,  as  can  be  seen  from  Fig.  6.14,  the  radius  of  maximum  actual 
wind  decreases  with  time.  Figs.  6.15  and  6.16  show  h(r)  and  £/ /.  As  noted  in  the  shallow- 
water  case  they  depict  the  adjustment  of  the  mass  field  and  the  wind  field.  The  adjustment 
here  is  in  the  same  sense.  As  the  potential  vorticity  field  evolves  the  contribution  from 
the  mass  field  increases  and  the  contribution  from  the  wind  field  decreases. 

The  results  from  the  upper  layer  show  some  interesting  features.  Fig.  6.17  shows 
v(r).  We  see  that  the  upper  level  anticyclone  intensifies  with  time  and  is  pushed  outward. 
This  is  also  evident  in  Fig.  6.18  which  shows  V(R).  However,  in  this  case  the  radius 
of  maximum  “wind”  decreases.  The  fluid  depth  and  vorticity  fields  (not  pictured)  show 
some  of  the  features  we  would  expect,  increasing  fluid  depth  and  decreasing  vorticity  with 
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Figure  6.13:  V  as  a  function  of  R  for  layer  one  corresponding  to  the  H(R,T)  fields  in  Fig. 

6.12. 
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Figure  6.16:  The  normalized  vorticity  (£//)  as  a  function  of  r  for  layer  one  corresponding 
to  the  H(R,T)  fields  in  Fig.  6.12. 
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Figure  6.17:  The  tangential  wind  (t>)  as  a  function  of  r  for  layer  two  corresponding  to 
T  =  1,2, . . . ,  7  days. 
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Figure  6.18:  Vasa  function  of  R  for  layer  two  corresponding  to  T  =  1,2, ...  ,7  days. 
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time  at  the  center.  However,  due  to  the  crude  lateral  boundary  condition  the  outer  region 
is  somewhat  anomalous.  We  get  increasing,  slightly  cyclonic  vorticity.  This  is  caused  by 
requiring  the  average  vorticity  be  zero. 


Chapter  7 


SUMMARY  AND  CONCLUSIONS 

We  have  attempted  to  solve  the  invertibility  principle  for  Ooyama’s  tropical  cyclone 
model  using  potential  vorticity  reasoning  and  potential  radius  coordinates.  Presentation 
of  potential  vorticity  and  potential  radius  in  a  shallow  water  model  showed  how  they  are 
used  to  an  advantage.  They  provide  stretching  in  regions  of  large  vorticity  for  better 
resolution  and  eliminate  the  radial  component  of  the  wind  from  the  problem. 

We  showed  that  the  invertibility  principle  can  be  written  and  solved  several  different 
ways.  Solving  for  $  or  h  using  the  shooting  method  works  very  well.  However,  these 
methods  can  not  be  generalized  to  two  layers,  because  the  secant  method  only  works  for 
problems  of  a  single  variable.  Solving  for  h  using  a  nonlinear  equation  solver  “black  box” 
also  worked  well.  Again  this  was  difficult  to  transport  to  the  full  model,  because  we  were 
not  able  to  uniquely  define  the  problem  in  this  case.  Solving  a  tridiagonal  matrix  equation 
for  V  looked  promising,  but  proved  disappointing.  This  method  produced  only  a  weak  vor¬ 
tex,  because  the  matrix  solver  was  too  sensitive  to  small  changes  in  the  vorticity.  Finally, 
solving  for  V  with  the  nonlinear  equation  solver  proved  the  best.  Although  somewhat 
slow,  this  produced  good  results  and  generalized  easily  to  two  layers. 

We  next  looked  at  Ooyama’s  model  and  how  it  solved  for  the  wind  field.  This  model 
is  based  on  gradient  balanced,  axisymmetric  flow  in  a  homogeneous  incompressible  fluid. 
Prediction  is  based  on  the  radial  mass  flux  and  involves  solving  a  system  of  three  cou¬ 
pled  pairs  of  equations.  The  results  show  a  vortex  which  evolves  into  a  tropical  cyclone 
after  approximately  six  days.  The  structure  and  evolution  are  in  strong  agreement  with 
observational  data. 
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Our  tropical  cyclone  model  incorporates  the  potential  radius  ideas  as  in  the  shallow 
water  case.  Although  we  include  friction  in  this  case,  the  advantages  of  coordinate  stretch¬ 
ing  and  elimination  of  the  radial  wind  component  remain.  The  invertibility  is  derived  and 
solved  using  the  nonlinear  equation  solver.  It  is  interesting  to  note  that  the  problem  has 
been  reduced  to  one  coupled  pair  of  equations. 

Through  this  analysis  we  have  shown  a  different  way  of  looking  at  tropical  cyclones. 
The  release  of  latent  heat  produces  positive  relative  vorticity  at  low  levels  and  negative 
relative  vorticity  at  upper  levels.  This  in  turn  produces  cyclonic  winds  at  low  levels  and 
anticyclonic  winds  at  upper  levels.  The  wind  and  mass  fields  mutually  adjust  to  the 
production  of  potential  vorticity  anomalies  by  latent  heat  release. 

This  work  has  shown  that  it  is  possible  to  create  a  tropical  cyclone  model  using 
potential  vorticity  as  the  predictive  variable.  However,  more  work  needs  to  be  done  to 
close  the  system  and  make  it  a  true  tropical  cyclone  model.  The  boundary  layer  must 
be  added  to  include  the  air-sea  interchange  of  moisture  and  surface  frictional  effects.  A 
parameterization  of  cumulus  convection  needs  to  be  added  to  more  accurately  describe 
the  release  of  latent  heat.  It  would  also  be  useful  to  make  the  model  less  restrictive  by 
allowing  the  user  to  specify  the  initial  wind  field,  the  sea  surface  temperature,  and  the 
latitude.  This  way  the  model  could  be  used  to  examine  different  initial  conditions  in 
different  locations. 
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